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i.  introduction 

A  Bpm,  smoothness  prion  approach  it  iatrodactd  in  this  paper  for  transfer  finctioa  esti¬ 
mation.  Jointly  stationary  input  nad  ontpat  data  it  nssumsd  to  hs  observed  in  the  presence  of 
additive  colored  noise.  The  method  it  pnrtkelariy  relevant  when  only  short  spaas  of  data  are 
available,  when  the  impalte  response  it  relatively  long- tailed  and  when  the  low  order  polynomial 
ARMAX  type  model  can  not  capture  the  tree  model  etrectnre.  The  method  is  illustrated  by  the 
analysis  of  the  Box-Jenkms  Series  J  data.  The  statistical  performance  of  the  method  it  explored  in 
Monte-Carlo  simulation  studies. 


The  models  in  Antrom  and  Bohlin  (1965)  and  Box  and  Jenkins  (1970),  are  the  classical 
parametric  time  domain  transfer  function  models.  In  that  method,  ARMAX  type  models  charac¬ 
terised  by  polynomial  operators  on  the  input,  the  output  and  the  observation  noise  are  fitted  to 
the  observed  input  and  output  data.  (The  observation  noise  in  the  Astrom  -Bohlin  model  is  MA 
noise.  It  is  A  A  noise  in  the  Box- Jenkins  model.)  That  method  requires  the  specification  of 
three  polynomial  operator  orders,  one  each  for  the  input,  output  and  noise  polynomials  and  the 
estimation  of  the  unknown  polynomials  coefficients  via  the  minimisation  of  a  performance  func¬ 
tional.  Typically  that  computation  is  achieved  by  a  computationally  costly  nonlinear  optimisation 
procedure.  In  such  procedures  it  is  only  feasible  to  search  for  solutions  over  low  polynomial  orders. 


Despite  the  fact  that  conventional  transfer  estimation  methods  have  been  extensively  used, 
the  influence  of  the  sampling  variability  in  the  polynomial  mode)  order  selection  on  the  transfer 
function  estimation  performance  remains  to  be  explored.  Another  objection  to  the  use  of  low 
order  polynomial  ARMAX  models  is  that  the  "parsimonious"  parametric  model  may  not  be  a  good 
characterisation  of  the  system  that  generated  the  data.  An  elaboration  of  this  objection  to  the 
conventional  parametric  modeling  method,  from  a  Bayesian  point  of  view,  is  that  the  conventional 
parametric  modeling  methods  can  not  yield  "correct"  models.  That  is,  if  there  is  information  in 
the  data  to  select,  by  some  best  model  order  selection  procedure,  an  ARMAX  p,q,r  model,  then 
there  is  also  information  in  the  data  to  select  alternative  ARMAX  p’.q'.r’  models.  There  the  best 
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Bagreeiao  mndsl  wyiim  that  the  transfer  hiKlion  ronpaud  with  dflfrrewt  mod*)  order*  be  ever- 
^  with  respect  to  the  likelihood  at  each  Itted  model  eod  the  prior  probabilities  at  the  model 
order*.  A  hoik*  (1979s)  it  perhaps  the  tot  example  of  a  Bayesian  time  series  modeliag  with 
averaging  of  the  compatatieaal  resalts  over  models  with  different  model  orders.  The  smoothness 
prises  method  of  transfer  faactioo  estimation  in  trod  seed  hers  obviates  the  9  stringent  parameter 
model  order  search  problem  in  conventional  transfer  estimation  procedures.  Oar  procedure  also 
uses  8  parameters.  One  is  a  1  dimensional  model  order  parameter,  the  other  two  are  smoothness 
differential  order  parameters.  We  demonstrate  that  the  values  of  those  parameters  are  not  critical 
in  determining  the  estimated  transfer  function  properties. 

The  Mayne-Firsoon  (1978)  three  stag*  least  squares  (SSLS)  procedure  was  developed  to 
avoid  the  costliness  of  the  maximisation  of  the  nonlinear  likelihood  for  the  transfer  function 
modeled  wirth  additive  M4  noise..  Almost  contemporaneous  alternative  linear  computational 
transfer  function  estimation  procedures  were  generalised  least  squares  (Clarke,  1987)  and  extended 
least  squares  (Paauska,  1968).  Durbin  (1961),  a  SSLS  procedure,  was  the  conceptual  predecessor 
of  the  Mayne-Firsoon  procedure.  A strom  and  Mayne  (1989)  and  Hannan  et  el.  (1986)  are  recur¬ 
sive  procedure  realisations  of  the  Mayne-Firsoon  SSLS  transfer  function  estimation  procedure. 
Other  recent  noteworthy  publications  on  or  related  to  transfer  function  estimation  include  Hannan 
and  Riasanen  (1989),  a  recursive  method  for  finding  model  orders,  and  Ljung  (1985).  a  study  of  the 
statistical  properties  of  time  and  frequency  domain  transfer  function  estimation  procedures. 

Jordinson  et  al.  (1970),  New  bold  (1979),  W  eg  man  (1980),  Jakemaa  and  Young  (1989),  and 
Kruc  et  al.  (1989)  are  examples  of  the  literature  on  statistical  regularisation  and  Bayesian 
smoothed  deconvolution  procedures  for  the  estimation  of  transfer  functions.  Applications  of  that 
literature  include  the  estimation  of  the  transfer  function  of  the  vascular  system,  applications  in 
radiology,  dispersive  relations  in  streams  etc.  This  activity  is  not  summarised  here  .  Our  own 
smoothness  priors  method  is  a  variation  on  that  Bayesian  theme. 
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b  m  method.  an  Mtk  order  impulse  response  between  input  end  output  plus  an  Mth  order 
uutoeegressive  (AR)  model  for  the  additive  noise  is  assumed,  with  M  "quite  large".  This  model  is 
equivalent  to  an  ARMAX  plus  white  noise  model.  We  assume  integrated  square  seroth  and  kth 
order  derivative  frequency  domain  smoothness  constraints  on  the  polynomial  operators  In  the 
least  squares  framework,  the  resultant  model  strikes  a  balance  between  the  infidelity  of  the  solution 
to  the  data  and  the  infidelity  of  the  solution  to  the  smoothness  constraints.  That  balance  or  tra¬ 
deoff  is  characterised  by  one  parameter  for  each  of  four  smoothness  constraints.  In  Bayesian  ter¬ 
minology,  those  are  referred  to  as  hyperparameters,  (Lindley  and  Smith,  1972).  The  likelihood  of 
the  hyperparameters  that  characterise  the  claw  of  smoothness  priors  is  maximised  to  yield  the 
best  transfer  function  model  with  the  best  data  dependent  priors. 

The  approach  taken  in  this  paper  is  an  application  of  our  frequency  domain  smoothness  pri¬ 
ors  AR  model  spectral  estimation  method,  Kitagawa  and  Cersch  (1985a).  Some  of  our  other 
time-domain  smoothness  priors  papers  on  the  modeling  of  nonstationary  time  series  are  Brotherton 
and  Gersch  (1981),  Kitagawa  (1981)  and  Kitagawa  and  Gersch  (1984,  1985b).  Additional  Bayesian 
methods  of  time  series  analysis  are  in  TIMSAC-84,  Akaike  et  el.  (1985).  Gersch  and  Kitagawa 
(1987),  is  a  review  of  our  smoothness  priors  modeling  of  time  series.  The  papers  by  Shiller  (1973) 
and  Akaike  (1979b,1980)  are  predecessors  to  our  own  work.  In  particular.  Akaike  (1980), 
motivated  our  interest  in  this  subject.  Additional  related  work,  that  is  better  known  in  the 
statistics  literature,  are  the  methods  of  regularisation,  Tikhonov  (1965),  and  the  maximised  penal¬ 
ised  likelihood  method  (I.J.  Good  1970,  and  Good  and  Gaskins  1980),  Wahba  (1977), (1982), (1983) 
and  O’Sullivan  (1987). 

The  smoothness  priors  method  of  transfer  function  estimation  is  treated  in  Section  2.  An 
example  of  the  smoothness  priors  analysis  of  the  Box-Jenkins  Series  J  data  is  shown  in  Section  3. 
Studies  of  the  statistical  performance  of  the  smoothness  priors  method  and  comparisons  of  the  SP 
and  SSLS  methods  of  transfer  estimation  are  also  in  Section  3.  An  interpretation  of  the  results, 
summary  and  discussion  in  Section  4  conclude  the  paper. 
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2  ANALYSIS 


la  Section  2.1,  several  features  of  a  Bayesian  model  for  linear  refression  are  shown.  A 
variety  of  other  analyses  of  the  Bayesian  model  are  shown  for  example  in  Zellner  (1971)  and 
Broemeling  (1985).  The  presentation  here  uses  assumptions  on  the  priors  of  the  parameter  vector 
similar  to  those  used  earlier  in  Lindley  and  Smith  (1972)  and  Akaike  (1980).  Illustrations  of 
smoothness  priors,  a  particularization  of  the  Bayesian  linear  model  analysis,  are  shown  in  applica¬ 
tions  to  the  time  series  analysis  problems  considered  by  Whittaker  (1923)  and  Shiller  (1973)  in 
Section  2.2.  Our  transfer  function  model  is  described  in  Section  2.3.  In  contrast  with  the  time 
domain  priors  on  the  model  coefficients  in  the  Whittaker  and  Shiller  problems,  in  Section  2.4  we 
show  an  application  of  smoothness  priors  in  the  frequency  domain,  in  transfer  function  estimation. 

2.1  A  BAYESIAN  LINEAR  REGRESSION  MODEL 

Consider  first,  the  linear  regression  model 

9  m  X0  +  t  .  (2.1.1) 

In  (2.1.1),  y  =(Vi>—iirjv)>  i»  *■  nxl  vector  of  observations,  9  is  a  pxl  parameter  vector,  X  is  an 
nxp  known  design  matrix  and  t  is  an  i.i.d.  nxl  random  vector  with  t  ~N(0,o*l).  Then,  the  con¬ 
ditional  data  distribution  is 

P(* I  XA*)  -  (2xo»)-/*exp| ~(y  -  X9  ) T[y  -  X6  )  j  .  (21.2) 

In  the  stochastic  regression  problem,  9  is  considered  to  be  a  random  vector  with  distribution  *(0). 
From  Bayes  theorem,  the  posterior  distribution  of  the  parameter  vector  9  is  proportional  to  the 
product  of  the  conditional  data  distribution  (the  likelihood),  p(y I  X,9,<r7),  and  the  prior  distribu¬ 
tion,  x(0). 

x(8|  y,a)  a  p(y |  X,9,o)*{9)  (2.1.3) 

Let  the  prior  distribution  of  9  be  D9  —  N(D90,Xil).  In  the  application  of  Bayesian  regression 
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methods  to  smoothness  prion  problem*,  we  shall  consider  the  special  case  of  #0  -  0  Also,  it  is  con¬ 
venient  to  introduce  the  parameterisation,  A  -  r/a.  la  that  case,  the  prior  distribution  on  3  is 


a/* 


*(#!  vr)  -  fDTD  exp{ --^tTDTDt 


,|-4< 

la3 


(2.1.4) 


r  is  referred  to  as  the  hyperparameter  of  the  prior  distribution,  (Lindley  and  Smith,  1972).  In  this 
conjugate  family  Bayesian  situation,  in  which  both  the  prion  and  conditional  data  distribution  are 
normally  distributed,  the  posterior  distribution  is  also  normally  distributed,  (Zellner  1971).  The 
mean  of  that  distribution  is  easily  computed  as  the  minimiser  of 

I  B 


(2.1.5) 


I 


If  r  were  known,  the  computational  problem  in  (2.1.5)  could  be  solved  by  an  ordinary  least 
squares  computation.  The  solution  for  i  is 

6  =  [xTX  +  t3DTD]  'xTy  (2.1  6) 

with  the  residual  sum  of  squares, 

SSE(r)  =  yTy  -  ST[xTX  +  i3Dtd]$.  (2.1.7) 

The  posterior  distribution  of  0,  r(6\  y ,r,«r)  is  a  proper  distribution,  therefore  the  likelihood 
for  the  unknown  parameter  r  can  be  determined  by 

Hr, a)  =  f  w[0\  y,r,o)i$  (2.1.8) 

*/  —00 


I.J.  Good  (1965)  referred  to  the  maximisation  of  (2.1.8)  as  a  Type  II  maximum  likelihood  method. 
Since  n(0\  y,r.<r)  is  normally  distributed,  (2.1.8)  can  be  expressed  in  the  closed  form,  (Akaike  1980), 

L(r, a)  =  (2x<rJ)-A,^|  i3DtD\  l'3\  XTX  +  r*DTD\  -,'Jexpj-^±.SS£(r)J  (2.1.9) 


The  maximum  likelihood  estimator  of  a1  is 
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d*  -  8SE{t)jN  . 


(2.1.10) 


h  is  convenient  to  work  with  -2  k>f  likelihood.  Using  (2.1.10)  in  (2.1.0)  yields 

(21.11) 

-2<ofL(rfd)  -  Mo$U  +  Mog(SS£(r)/y)  +  log)  XTX  +  i*DTD\  -  log!  r *DrD  +  N. 

This  is  the  bask  relation  that  we  use  in  oar  smoothness  priors  least  squares  analysis.  A  practical 
way  to  determine  the  value  of  f3  for  wkkh  the  -2iog-likelihood  is  minimised,  is  to  compute  the 
likelihood  for  discrete  values  of  r3  and  search  the  discrete  ~21og  likelihood-hyperparameter  space  for 
the  minimum.  If  there  are  more  than  say  2  hyperparameters,  it  might  be  more  expeditious  to  use 
a  gradient  search  algorithm  to  determine  the  hyperparameters  that  maximise  the  likelihood. 
Akaike  (1980)  demonstrated  the  first  practical  use  of  the  likelihood  of  the  Bayesian  model  and  the 
use  of  the  likelihood  of  the  hyper-parameters,  as  a  measure  of  the  goodness  of  fit  of  a  model  to  data. 

Several  other  facets  of  stochastic  regression  may  be  of  interest.  The  solution  of  the  ordinary 
least  squares  regression  problem  in  (2.1.1)  is 

hs  *  (XTX )-«Xry  (2.1.12) 

Matrix  algebra  yields 

9  -  (XTX  +  i*DTD)-l(XrXiLS  -r  t*DTD09).  (2.1. IS) 

That  is,  the  posterior  parameter  estimate  is  a  weighted  sum  of  the  least  squares  solution  and  the 
prior  mean,  09.  Let  0,  be  the  true  value  of  the  parameter  vector  0.  Then,  direct  evaluation  of  the 
mean  square  parameter  vector  error,  MSE{i)  «  Vor(^)  +  E{0—  0,)T  E((&-9t),  yields  the  result 

MSE[9)^MSE{9u)  (2.1.14) 

iff 

tr{XTX)-'2tr{XT X  +  i*DTD)-'  +  rJ(fi,-fio)r0r(*7'*  *  t*Dr  D)-*D(9,-0O)  . 

The  first  term  on  the  RHS  of  (2.1.14)  is  not  larger  than  the  LHS  of  (2.1. IS).  Depending  upon  how 
close  09  is  to  0„  the  MSE{i)  may  or  may  not  be  less  than  MSE[0ls).  The  Bayesian  method 
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minimises  expected  lorn.  Therefore  the  expected  value  of  MSE{9)  will  be  less  than  or  equal  to  the 
expected  value  of  MSE(its). 


1.9  SOME  EXAMPLES  OP  SMOOTHNESS  PRIORS  MODELING 

Two  of  the  earliest  smoothness  priors  problems  are  illustrated  here.  We  refer  to  those  as  the 
Whittaker  problem,  (Whittaker  1923),  and  the  Shiller  problem  (Shiller  1973), 

The  Whittaker  Problem:  in  the  problem  treated  by  Whittaker;  the  observations  ,,.V 

are  given.  They  are  assumed  to  consist  of  the  sum  of  a  "smooth"  function  and  observation  noise, 

(22.1) 

The  problem  is  to  estimate  the  unknown  /,,n  =  l,...,Af.  In  a  time  series  interpretation  of  this  prob¬ 
lem,  is  the  trend  of  a  nonstationary  mean  time  series.  A  typical  approach  to  this 

problem  is  to  use  a  class  of  parametric  models.  The  quality  of  the  analysis  is  completely  depen¬ 
dent  upon  the  appropriateness  of  the  assumed  model  class.  A  flexible  model  is  desirable.  In  this 
context,  Whittaker  suggested  that  the  solution  balance  a  tradeoff  of  goodness  of  fit  to  the  data  and 

goodness  of  fit  to  a  smoothness  criterion.  This  idea  was  realised  by  determining  the  /„,n  =  l . JV 

to  minimise 

IE  iln-fn)*  +  **E  (V7.)*l  (2.2.2) 

for  some  appropriately  chosen  smoothness  tradeoff  parameter  p*.  In  (2.2.2)  y*/»  expresses  a  kth- 
order  difference  constraint  on  the  solution  f,  with  y/i,  =  Jn  -  /,_ |.  y/*  =  V(V/«),  etc. 
(Whittaker’s  original  solution  was  not  expressed  in  a  Bayesian  context.  Whittaker  and  Robinson. 
1924  does  invoke  a  Bayesian  interpretation  of  this  problem.) 

The  properties  of  the  solution  to  the  problem  in  (2.2. 1  )-(2  2.2)  are  clear.  If  ^J  =  0,  /„  =  y, 
and  the  solution  is  a  replica  of  the  observations.  As  ft1  becomes  increasingly  large,  the  smoothness 
constraint  dominates  the  solution  and  the  solution  satisfies  a  kth  order  constrain!.  For  large  /t3 
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Mid  fc-1,  the  solution  is  a  constant,  for  fc- 2,  it  is  a  straight  line  etc..  Whittaker  left  the  choice  of 
ft*  to  the  invesi  i gator 


Prom  the  Bayesian  point  of  view,  the  difference  equation  constraints  on  the  parameter  vector 
problem  are  stochastic.  That  is,  “  w».  w>th  wn  assumed  to  be  an  i.i.d.  normally  distributed 
sero-mean  sequence  with  unknown  variance  r2  .  For  example  for  i”  1  and  k- 2  those  constraints 


/.  *  /.-i  + 

/•  “  2/«- 1  -  /•-:  + 


12.2.3) 


Corresponding  to  the  matrix  D  in  (2.1.8),  for  difference  orders  4*1  and  i  =  2  respectively,  the 
smoothness  constraints  can  be  expressed  in  terms  of  the  NxN  constraint  matrices  Z?,  and  D2, 


D,- 


-1  1 
-1  I 


-1  1 


-9  9 
1  -2  1 


0 

1  -2  1 


(2.2.4) 


In  (2.2.4)  a  and  9  are  small  numbers  that  are  chosen  to  satisfy  initial  conditions.  (An  alternative 
to  the  ad  hockery  in  specifying  a  and  9  in  (2.2.4)  is  to  estimate  f0  and  /,  by  a  maximum  likeli¬ 
hood  method.) 


We  use  the  parameterization  n=o,r.  Therefore.  has  a  noise-to-signal  interpretation. 
Larger  r  corresponds  to  smoother  trends.  For  fixed  k  and  fixed  i2  the  least  squares  solution  can  be 
expressed  in  the  form  of  (2.1.5).  The  matrix  X  and  the  parameter  vector  9  in  (2.1.5)  are  replaced 
by  the  identity  matrix  /  and  the  parameter  vector  /  —  (/,.. ../*)•  Then  for  example,  with  k  =  2  and 
D  •  Dj.  the  solution  {/..n  «  1,...,;V)  satisfies 


From  (2.1.6), the  solution  to  (2.2.5),  with  D=Di,  is 

/  =  (/  +  tDtD^'y  , 

and  the  value  of  SSE(r)  is  given  by  (2.1.8)  with  6=f,X=I,D*sDi. 
likelihood  for  this  problem  is: 


(2.2.5) 


(2.2.6) 

The  minimized  value  of  -21og 

(2.2.7) 


-2logL(f,a)  =  Nloglx  +  Nlog(-^SSE(r))  -  tog\  ?D?Dt  +  /|  -  log\  +  N. 


The  Shiller  Problem:  The  problem  treated  by  Shiller  is  the  estimation  of  the  distributed  lag 
(impulse  response),  given  the  jointly  stationary  time  series  observations,  {yn,xn;  n  =  The 

distributed  lag  model  is, 

U 

Vn  m  2  d"  *»  (2.2.8) 

m*0 

Frequently  and  also  in  the  case  of  the  data  analyzed  by  Shiller,  econometric  data  is  short  duration. 
As  a  result,  econometricians  have  been  motivated  to  Bayesian  analyses.  They  assume  a  prior  dis¬ 
tribution  on  the  parameters  of  the  model  and  thereby  increase  the  effective  data  length.  The 
smoothness  priors  assumed  on  the  distributed  lag  coefficients  by  Shiller  were  of  the  form, 
y*hm  =*wm,  with  {u>m  ,m=0,...,Af  },  a  zero  mean,  normally  distributed  zero  mean  i.i.d  sequence. 
Those  are  the  same  priors  assumed  for  the  Whitaker  example.  We  take  h0  to  be  0  and  for  simpli¬ 
city  consider  the  initial  conditions,  as  known.  Then,  the  computational  problem  for  the 

smoothness  priors  distributed  lag  parameters,  is  as  in  (2.1.6).  In  the  application  of  (2.1.6)  to  the 
Shiller  problem,  the  parameter  vector  9=h,  and  the  Nxm  matrix  X  as  given  below  in  (2.2.9)  and 
the  matrix  D  either  Dt  or  £)2  45  *n  (2.2.4). 


9 


IK-. 


p 


*0  •  •  •  X\-H 

9i 

*'•  • 

93 

Aj 

jr  - 

• 

.  9  - 

.h  « 

• 

*M-\  •  ■  •  *H-ii 

9h 

A* 

(*•*•») 


3.3  A  TRANSFER  FUNCTION  MODEL 

Assume  that  inpat /oat pat  jointly  stationary  time  series  data  za,pn  n-1  is  observed. 
Assume  that  the  output  ya  is  observed  in  the  presence  of  additive  colored  noise,  «ra.  Consider  a 
representation  of  the  input/outpat  plus  noise  in  the  impulse  response  plus  colored  noise  form. 


■>•1 


-  E  -*■  •« 

•■I 


(2.3.1a) 


(2.3.1b) 


For  convenience  in  (2.3.1b),  a.  is  assumed  to  be  a  Gaussian  sero-meaa  nncorrelated  sequence  with 
unknown  variance  rj.  In  (2.3.1a)  bm  is  an  impulse  response  sequence  and  w,  is  assumed  to  be  in 
AR  model  form. 


Using  the  assumed  stationarity,  (2.3.1a)  yields 


*«-/  "  !«-»  E ■ 

i-i 


Substituting  the  expression  for  wn_;  into  (2.3.1k)  yields  the  model 

00  00 

*  » 1  m  « I 

with 

em  *  «». 

*»-l 


(2.3.2) 


(2.3.3) 


(2.3.4) 


*m  “  -  E  m“1- 
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Eqmtioa  (2.3.3)  is  an  ARMAX  model  with  additive  white  noise  un.  The  models  in 
(2.3.1a), (2.3.1b)  are  estimated  using  (2.3.3)  and  (2.3.4).  The  infinite  order  transfer  function  model 
in  (2.3.3)  is  approximated  by  finite  transfer  function  model 

w  M 

f»  “  S  c«k«~«  ^  *n-m  +  a*  i  (2.3.5) 

■•I  «*»t 

with  M  assumed  to  be  "large".  (The  choice  of  M  may  be  determined  by  the  maximization  of  a 
likelihood  and  Akaike’s  A1C.)  The  coefficients  em,dm  m  =  l are  directly  estimated  by  the 
Bayesian  procedure  described  in  the  following  section.  The  estimates  of  the  coefficients  of  the 
model  (2.3.1a)  are  then  obtained  by  the  formulas 

°ib  “  w*“l,...)A/  (2.3.6) 

am  *  0,  m=M+l,... 


W-l 

6m  *  dm  +  £)  flj6m— m=l,...,M 

i-i 

M 

bm  =  Bi-Af+1,M+2|...  . 

•»i 


From  (2.3.1a),  the  frequency  response  function  from  the  input  z„  to  the  output  y„  can  be 
obtained  from 


*(/)  “  £  &m«XPi-2*t>»] 


where  «*— — 1.  The  power  spectrum  of  the  noise  wn  is  given  by 


5(/)  -  - - - ^ - 

I  1  -  £  «»exP!_2x»/m]| 


where  a *  is  the  innovations  variance  of  the  estimated  model  (2.3.5). 
Identify  the  quantities  C(f )  and  D[J ) 


(2.3.7) 


(2.3.8) 
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(2.5.9) 


introduce  the  ttro-tk  derivative  smoothness  constraints 

*0  -  C  Cm  Ui  =  1  +  S  4  (2.4. S) 

Go-  f_j  ix/MVr-js « 

Let  the  differential  orders  for  the  numerator  and  denominator  polynomials  to  be  Jk,,  and  k3  respec¬ 
tively.  With  these  "frequency  domain"  priors  we  then  have  the  constrained  least  squares  problem 
which  for  fixed  values  of  tl,ki  and  rf,  determines  the  {en,dm,  m  =  l,...,M)  that  minimises 


(2-4.4) 


N 

E(v» 


n-1 


E  ~  E  tmtn-m?  +  E  I»l4  +  ^^‘‘4! 

**1  «■!  m»l 


E  [rfd*  + 


In  (2.4.4)  rf  /=1,....4  are  the  tradeoff  parameters.  By  a  proper  choice  of  the  tradeoff  parameters, 
our  estimate  of  the  model  parameters,  {cm,dml  m  =  l,...,M),  balance  a  tradeoff  between  infidelity  of 
the  transfer  function  solution  to  the  data  and  infidelity  to  the  smoothness  constraints. 


2.5  THE  SMOOTHNESS  PRIORS  LEAST  SQUARES  PROBLEM,  DETERMINING 
THE  TRADEOFF  PARAMETERS 


As  indicated  in  Section  2.1,  the  constrained  least  squares  problem  has  a  Bayesian  interpreta¬ 
tion  which  facilitates  the  determination  of  the  tradeoff  parameters  in  the  criterion.  We  apply  that 
approach  to  the  particular  smoothness  prior  transfer  function  problem  at  hand. 

In  detail,  the  minimisation  of  (2.4.4)  is  equivalent  to  the  maximisation  of 


(2.5.1) 


exp< 


ttE1v«  -  E e mV ■- 

2<*  »•]  m- 1 


M 

-  E 

m*  1 


.-J*|«xpj“ ■  E  (*f  +  fs»n2*,)e*|exp|-^  £  (rf  +  rfm“: 


In  that  form,  the  constrained  least  squares  solution  has*  a  Bayesian  interpretation  as  the  maximum 
a  posteriori  estimate  of  the  model  with  the  data  distribution 


)4 


is 


;*  ’  • 

J-.-*  *:., 


P' 


I": 


J»(fl 

ud  the  prior  distribution 


XA*)  -  (2^rw/J«p|-jEi».  -  E  «-fc-  -  E  O.-J’l  (2.5.2) 

\  2<r  «*)  m«]  m»  1  J 


ir(*|  r,«r)  -  (2»«t*)-"|  DTD\  ^cxp\^\«T D1  D9\ 

2<r 


(2.5.3) 


where  r  denote*  the  vector  of  hyperp ammeters  r,,...r4)  and  6  denotes  the  model  parameters, 
. M.  In  (2.5.3) 


W-Mf)1'1 

(ri+ri)'/* 

Cl 

*x 

(r?+2,4‘r})>/> 

. 

eH 

du 

(2.5.4) 


From  (2.5.2)  and  (2.5.3)  it  follows  that 


J»(fl  X,$,o)w(9\  r,o) 

-I 


(2.5.5) 


(2«r*)-<"+MW*|  DTD\  ‘/»expi— i-(#-#)r(Xrjr  +  DTD)-'[i  -  #))**p(^±SS(r)], 

2 r  2  o3 


where 


Vo 

*0 

■  ■  9t-u 

Z|tf 

9i 

ft 

■  ■  9l-U 

*»-* 

Hi 

• 

,  *  - 

9h-i 

*N-1 

■  9N-U 

*N-M 

9h 

(2.5.6) 


and 


i  -  (XTX  +  DtD)  'Xti  , 
S5(r)  -  *r*  -  iT{ XTX  +  DTD)i 


(2.5.7) 


In  (2.5.6),  the  initial  condition  data  {*„>,,  t  «  -M,-A/+1,...,0)  are  assumed  known.  Integrating 
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(3,5.5)  with  respect  to  the  parameter  vector  yields  the  likelihood  of  the  hyperparameters, 


L(r,o)  -  (1*<t)-"l3\  DTD\  ,/s|  XTX  *  DTD  l/aexp!  ■  55(rJ 

la1 


(2.5.S) 


Then 


logI(r,o)  =  -^log2»<r*  +  -i-lo*|  DrD\  -^-logl  XTX  +  DTD\  -  -^j-SS(r),  (2.5.9) 


with  the  maximum  likelihood  estimate  of  a3  given  by 


(2.5.10) 


Substituting  the  estimated  value  of  a3  from  (2.5.10)  into  (2.5.9)  yields 

logI(r,o)  =  -^-log2rora  +  ylog|  DTD\  -  ylog|  XT X  +  DrD\  -  y  ,  (2.5.11) 

which  is  to  be  maximised  to  obtain  the  maximum  likelihood  estimates  of  tJ  ;»1,...,4.  The  likeli¬ 
hood  for  the  hyperparameters  is  maximised  via  a  Davidon-Fletcher-Powell  gradient  search  algo¬ 
rithm.  That  algorithm  is  exterior  to  a  Householder  transformation  least  squares  solution  of  the 
constrained  least  squares  problem. 

We  use  the  AIC  statistic,  (Akaike  1973), 

AIC  -  -2logL(r,o)  +  2(numher  of  parameters  estimated)  (2.5.12) 

•  miofiwa3  +  l)-k>g|  DTD\  +  log |  XTX  +  DTD\  +  2(numker  of  parameters  estimated), 

to  determine  the  order  M  for  the  transfer  function.  Akaike  (1980)  referred  to  (2.5.12)  as  the  likeli¬ 
hood  of  the  Bayesian  model.  The  analysis  indicated  here  is  referred  to  as  a  "quasi- Bayesian  " 
analaysis.  A  more  thorough  Bayesian  solution  of  the  transfer  function  estimation  problem,  would 
require  that  priors  be  specifled  on  the  model  order  M.  A  completely  orthodox  Bayesian  analysis 
would  require  priors  on  the  model  order  and  parameters  of  the  stochastic  input  to  the  system. 


S.  ANALYSIS  OP  BOX-JENKINS  SERIES  J  DATA  AND  MONTE  CARLO 
RESULTS 

la  this  section  the  results  of  the  transfer  function  analysis  of  the  Box-Jenkins  Series  J  gas  fur¬ 
nace  data  by  the  smoothness  prion  and  Box-Jenkins  methods  are  described  and  compared.  Some 
properties  of  the  smoothness  prion  method  of  analysis  are  shown.  Also,  we  show  the  results  of 
Monte  Carlo  studies  of  the  statistical  performance  of  the  smoothness  prion,  (SP),  method  and  the 
SSLS  asymptotically  maximum  likelihood  method  of  Mayne  and  Firsoon  ,(1978),  based  on  models 
derived  from  the  BJ  series  J  data. 

The  input  output  BJ  series  J  data  are  shown  in  Figure  1.  The  variances  of  the  input  and 
output  data  are  1.14727  and  10.25357  respectively.  For  illustrative  purposes,  the  data  shown  in 
Figure  1  was  normalised  to  have  sero  mean  and  the  same  variance.  (Inverting  the  output  data 
and  superimposing  it  over  the  input  data  reveal  the  output  to  be  a  delayed-low  pass  filtered  ver¬ 
sion  of  the  input.)  The  generic  Box-Jenkins  transfer  function  model  is 

ft.  -  «ift.-i  +  •  •  ■  -  «,ft.~,  +  M.-1-rf  -  +  V„-f -4  +  (3.1) 

*»  •  +  •  •  +  +  «„• 

In  (S.l)  (*,,»., w,,  ,n-l,...,JV)  are  respectively  the  observed  input  and  output  and  the  unobserved 
added  noise.  Also  in  (S.l),  {e„  u-1  ,...,/V)  is  a  norma)  i"ro  mean  i.i.d.  random  variable  with 
variance  o*.  The  dimensional  parameters  of  the  BJ  model  are  <f= 2,  p= 2,  g=3,  r«2.  The  pub¬ 
lished  vectors  of  BJ  model  coefficients  are:  «  =  (0.57,0.01);i  =  (0.53, -0.37, -0.51);e-(l. 53, -0.63), 
o*  -  0.05058,  (Box-Jenkins,  Section  11.4).  For  the  A1C  optimal,  *,=4,kj=2,  SP  model,  the 
dimensional  parameters  are  p=g=r= 4.  The  d=0  model  is  the  AIC  best  shift  parameter  model.  (In 
this  data  example,  the  higher  order  SP  model  automatically  accounts  for  the  delay  between  input 
and  output  data,  without  requiring  an  additional  non-sero  d  parameter.)  The  SP  a. 6  polynomial 
coefficients  are  s-(1.58824. -0.70509.-0.13198,0.14861),  6  =  (0.17090, -0  43813,-0.17497,0.08608). 
The  vector  of  smoothness  priors  tradeoff  parameters  was  ^=(0.00886,0.10305,0.71914.0.10686). 
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Ttblc  1  shows  the  values  of  the  AIC  for  the  differential  orders  k,,l2=l,...,4  for  the  order  M= 4 
node). 

INSERT  FIGURE  1  HERE:  BOX-JENKINS  SERIES  J  DATA 


Table  1.  AIC’ 

s  d=0,  M=4,  S 

P  Model,  Part 

imetric  in  kl  a 

nd  k2. 

k2 

kl  =  1 

kl  *  2 

kl  =  3 

kl  =  4 

kl=5 

1 

46.886 

45.697 

45.663 

46.345 

47.375 

2 

46.486 

46.186 

45.316 

45.995 

46.630 

3 

46.273 

45.031 

45.120 

45.794 

46.436 

4 

46.672 

45.013 

45.096 

45.769 

46.414 

5 

46.410 

45.072 

45.072 

45.721 

46.388 

In  Figure  2,  the  original  Box-Jenkins  Series  J  output  data,  and  the  SP  model  tracked  output 
data  are  shown  superimposed.  A  vertical  scale  displaced  version  of  the  difference  between  the  ori¬ 
ginal  and  tracked  data  is  also  shown  in  Figure  2.  The  tracked  output  data  is  computed  by  passing 
the  input  through  the  estimated  model.  The  appearance  of  the  S =  294  tracked  data  version  of  the 
BJ  modeled  data,  incorporating  the  d-2  parameter,  appears  very  similar  to  that  of  the  SP 
modeled  data  and  is  not  shown  here. 

In  the  sense  of  minimum  mean  square  tracking  error,  the  performance  of  the  AIC  optimum 
SP  and  BJ  models  were  similar.  The  variances  of  the  tracking  error  for  the  BJ  and  SP  models  (the 
sums  of  squares  of  the  residuals,  SSE),  were  0.70187  and  0.68662  respectively.  The  ratio  of  the 
relative  variance  of  the  tracking  error  to  the  variance  of  the  true  output  was  0.06798  and  0.06696 
for  the  BJ  and  SP  model  respectively. 


INSERT  FIGURE  2  HERE:  TRACKED  DATA 
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The  impulse  response,  transfer  function  amplitude  response,  phase  response  and  power  spec* 
tra)  densities  (psds)  of  the  residuals  associated  with  the  BJ  and  SP  models  are  shown  in  Figure  S. 
(Compare  the  transfer  function  and  residual  spectrum  of  a  windowed  periodogram  analysis  of  the 
B-J  data,  Jenkins  and  Watts,  p446.)  The  residual  psds  were  computed  from  Householder 
transformation*  Akaike  AIC  criterion  AR  models.  The  coefficients  of  the  corresponding  AIC  best 
ARt  models  of  the  BJ  and  SP  residuals  were  {1.53458,-0.55879,-0.21378,0.16280}  and 
{1.58329,-0.65524,-0.17195,0.17461}  respectively  with  corresponding  innovations  variances 
0.05995  and  0.05755.  In  Figure  3,  the  SP  and  BJ  modeled  psds  of  the  residuals  are  almost  identi¬ 
cal.  The  AR-A1C  model  of  the  residual  of  the  SP  modeled  BJ  Series  J  data  shown  in  Figure  3  is 
quite  similar  in  appearance  to  that  obtained  automatically  by  the  SP  modeling  procedure  and 
computed  directly  by  equation  (2.3.8). 

Also  in  Figure  3,  after  the  first  3  time  points,  when  the  impulse  response  for  the  BJ  model  is 
sero,  the  impulse  responses  of  the  BJ  model  and  SP  model  appear  similar.  The  SP  model  impulse 
response  is  smoother  than  the  BJ  impulse  response.  The  initial  non  sero  going  part  of  the  SP 
impulse  response  is  a  consequence  of  the  fact  that  the  optimal  SP  model  delay  parameter  d  is  sero. 

The  BJ  modeled  transfer  function  and  phase  function  versus  frequency  each  have  some  rela¬ 
tively  abrupt  kinks  in  their  responses  as  compared  to  those  for  the  SP  modeled  results. 

INSERT  FIGURE  3  HERE:  IMPULSE  RESPONSE,  TRANSFER  FCN  Sc  PHASE  Sc 
NOISE  PSD’s  SP  A  BJ  MODELS 

Some  comments  on  the  stationarity  of  the  Box-Jenkins  Series  J  data  are  in  order  here.  In 
Figure  2,  the  true  BJ  output  data  and  the  SP  modeled  data  appear  more  discrepant  in  the  latter 
part  than  the  earlier  part  of  the  time  series.  Also  in  Figure  2,  there  are  relatively  large  excursions 
in  the  latter  part  of  the  residual  time  series.  That  evidence  suggests  that  the  Series  J  data  is  nons¬ 
tationary.  To  examine  that  conjecture,  we  examined  the  residuals  of  the  AR  modeled  tracked 
data,  fitted  SP  models  to  the  first  200  data  points  and  to  the  first  75  data  points  of  the  Series  J 
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data  ud  examined  the  tracking  behavior  of  those  models.  The  properties  of  the  BJ  series  J  data 
do  appear  to  change  slightly  after  n  =  225. 

The  impulse  response,  noise  spectrum,  transfer  function  and  phase  function  for  the 
fci“l,tj*l  SP  models  for  n  =200  and  n =75  data  are  shown  superimposed  in  Figure  4A  and  4B 
respectively.  As  expected  from  Figure  2,  the  properties  of  the  SP  modeled  n=296  and  n=200  series 
J  data  are  quite  similar.  The  properties  of  the  SP  modeled  n  =  296  and  n  =  75  data  are  also  quite 
similar.  The  conjecture  that  the  SP  modeling  method  might  be  reasonable  for  relatively  short 
length  data  spans  is  supported  by  the  evidence  shown  in  Figure  4.  The  apparent  property  of  the 
SP  procedure  to  yield  reliable  parameter  estimates  with  relatively  short  data  length  time  series  is  a 
consequence  of  the  assumption  of  priors  on  the  model  parameters.  The  priors  are  equivalent  to  the 
observation  of  additional  data. 

For  completeness,  the  a,b  polynomial  coefficients  corresponding  to  SP  Af=4,it,  =  1,^=1 
modeled  data  are  respectively:  {1.13620,-0.15871,-0.29041,0.13128}, 

{1.08946,-0.24032,-0.55796.0.11972}  for  the  n=200  data  and 

{1.25647,-0.20495.-0.35709.0.16760},  {1.18340,-0.54000,-0.36671.0.28644}  for  the  n  =  75  data. 
The  respective  residual  variances  of  the  SP  modeled  Af“4,  n=200  and  n  =  75  data  point  models 
were  0.09900  and  0.04629  respectively.  The  corresponding  relative  tracking  variance  ratios  were 
0.01078  and  0.00620. 

INSERT  FIGURE  4  HERE:  IMPULSE  RESPONSE  etc  n=200  and  n=75  Models 

The  effects  of  the  choice  of  model  order.  Af,  and  the  differential  orders  , A3  on  the  impulse 
response,  the  noise  spectrum  and  the  transfer  function  amplitude  and  phase  of  the  SP  model  are 
also  of  interest. 

First,  to  illustrate  the  effect  of  model  order  Af  on  transfer  function  model  properties,  graphs 
of  impulse  response,  amplitude  and  phase  response  are  shown  superimposed  in  Figure  5A,B,C  for 
the  optimal  SP  A/=4  model  and  the  likelihood  best  SP  models  of  orders  A/=  10.20  and  Af=30.  The 


graphical  results  for  the  higher  order  SP  models  wiggle  only  slightly  around  those  for  the  order 
M=4  model.  Those  results  indicate  that,  on  the  provision  that  the  model  order  is  sufficiently 
large,  the  specification  of  the  order  of  the  SP  model  does  not  very  critically  influence  the  transfer 
function  characteristics.  For  completeness,  some  of  the  computed  results  for  those  models  are: 
M-10:  *1-1,*,=4,A/C=47.789;  M-20:  =4,^4^-50.224;  M-SO:  *,  =  2, *,=4, ,4/056.886. 

The  similarity  in  the  appearance  of  the  model  properties  for  the  different  model  orders  is  compati¬ 
ble  with  the  similarity  in  the  values  of  the  AIC  for  the  different  models. 

The  values  of  the  M=  10  optimal  SP  model,  a  and  b  polynomial  coefficients  are: 
a  «  {1.56980,-0.70201,-0.02424,0.06248.0.01102,-0.00000,-0.00152,-0.00046,0.00004,0.0010}, 
b  -  {0.15110,-0.38420,-0.17713,-0.00672.0.05029,0.08646,0.04176.0.00577,-0.01716,-0.00717}. 
The  pattern  of  a, 6  polynomial  coefficients  is  similar  for  the  larger  order  SP  models.  The  tapering 
toward  sero  values  effect  of  the  smoothness  priors  constraints  on  the  model  parameters,  particu¬ 
larly  on  the  higher  order  a  polynomial  parameters  and  the  relatively  long  tail  b  model  parameters 
helps  explain  the  similarity  of  the  M-4, 10,20  and  M-30  model  properties.  The  b  model  parame¬ 
ters  in  the  numerator  of  the  rational  polynomial  description  of  the  model  do  not  have  as  dramatic 
an  effect  as  do  the  a  polynomial  denominator  polynomial  parameters  on  the  model  properties. 
The  long  tail  b  polynomial  parameters  and  the  short  tail  a  polynomial  parameters  are  well  approx¬ 
imated  by  the  SP  M  =  4  model. 

For  the  purposes  of  comparison,  we  also  fitted  ordinary  least  squares,  (OLS),  models  of  ord¬ 
ers  Af-5,10  and  M-20  to  the  BJ  series  J  data.  Graphical  results  of  the  impulse  response,  transfer 
function  and  phase  function  for  those  models  are  shown  in  Figures  5D,E,F.  The  M=  5  LS  model 
properties  are  very  similar  to  the  optimal  SP  model  properties.  (The  relative  variance  of  the  track¬ 
ing  error  was  0.06706.  The  OLS  M=5  model  is  actually  a  superior  model  of  the  BJ  series  J  data 
than  the  original  BJ  model.)  The  computed  properties  of  the  OLS  M=  10  and  M=  20  models  wig¬ 
gle  a  lot  more  around  the  SP  M-4  model  than  the  SP  M=  30  model.  This  is  very  clear  evidence 
that  the  SP  model  properties  are  relatively  insensitive  to  model  order  in  comparison  with  other, 
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INSERT  FIGURE  f  HERE:  EFFECT  OF  MODEL:  ORDER  M 
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li**’  «ahaff  (*,4.i«l-lfJ-M  6)  ware  aaeasaad  known  la  modeling  data  with  aa  SP  model 

ordar  If  see  na*amh  Mhe  |l,.g.i>l.  If)  as  ataatl  roadatoa*  aad  model  the  data  on  the 
l**®***f  N -  If  data  pomea  The  hhahhaod  a  actually  competed  for  the  last  .V-|f  observation*, 
(hf-i-  1*1  kn  the*  cane,  modsh  *f  diSmat  If  order*  are  modeled  os  different  data  aad  it  is  not 
appropriate  to  nee  the  AJC  to  dataphil  beta  sea  me  del*  A  formula  that  permit*  the  AIC’»  of 
model*  of  ddlwai  order*  to  hr  rrwiparif  a. 


-  3( another  e/  parameter*  estimated).  (3  2) 

The  formala  in  (3.3)  ie  reasonable  seder  the  assamptioe  that  the  data  is  stationary,  i.e.  the  pro¬ 
perties  of  the  Irst  M  vales  of  the  data  do  not  change  very  much  with  different  values  of  M 


31 


Monte  Carlo  Results 


Monte  Carlo  simulation  studies  were  performed  to  explore  the  statistical  performance  of  the 
SP  method  of  transfer  function  estimation.  Some  comparisons  of  the  SP  and  Mayne  and  Firzoon 
(1978)  SSLS  procedure  were  also  done.  We  chose  to  compare  the  SP  with  the  SSLS  procedure 
because  the  latter  procedure  has  asymptotic  MLE  properties  and  is  easy  to  program.  The  principal 
topics  of  interest  are  the  bias  and  MSE  parameter  estimation  properties  and  transfer  function 
confidence  interval  properties  of  the  SP  and  the  SSLS  procedures.  We  show  computational  results 
of  simulation  studies  of  the  SP  and  SSLS  models  with  additive  AR  and  additive  MA  observation 
noises. 

Consider  the  transfer  function  model 

yn  =  *  •  •  +  afy<i-p  +  Mn-l  +  •  •  •  +  hiXn-i  +  (3-3) 

f 

with  vn  an  added  noise  process.  The  MAr  noise  model  is:  t„  =  £  e>nun~m  where  {u„}  is  a  zero 

m  =0 

f 

mean  i.i.d.  process.  The  ARr  noise  model  is  v„  =  V  cmt»n_m  *  un  .  where  again  {u„}  is  a  zero 

m  =  1 

mean  i.i.d.  process.  The  AR  observation  noise  model  is  the  Box-Jenkins  model.  The  MA  observa¬ 
tion  noise  model  was  used  in  Astrom  and  Bohlin  (1965).  Since  then  it  has  been  used  extensively  in 
engineering  applications. 

The  SSLS  procedure  was  developed  as  an  alternative  to  the  comptationally  extensive  max¬ 
imum  likelihood  method  for  the  MA  observation  noise  model.  For  convenience,  the  SSLS  pro¬ 
cedure  is  as  follows: 

Let  a,b.c  denote  the  AR,  MA  and  added  MA  noise  polynomials  respectively  in  equation  (3.3). 

i) lsing  least  squares  (LS),  fit  a  "long"  a.b  polynomial  model,  to  the  {z„,yn,  n  =  l,...,A'}  input- 
output  data  and  compute  the  residual  time  series. 

ii)  Fix  the  orders  of  the  a.b  and  c  polynomials  to  their  final  model  orders  and  use  the  original 
input-output  data  and  the  residuals  from  stage  i)  to  estimate  the  a.b  and  c  polynomial  coefficients 


by  LS. 

iii)  Prewhiten  the  input  and  output  data  using  the  inverse  of  the  c  polynomial  determined  in  stage 
ii)  and  estimate  the  fixed  order  a,b  polynomials  coefficients  by  LS. 

We  fit  the  SSLS  model  to  the  original  Box-Jenkins  Series  J  data  in  order  to  verify  the 
relevance  of  that  procedure  for  a  comparison  of  results  AR  observation  noise  Monte  Carlo  study. 
An  first  stage  SSLS  a,b  polynomial  model  order  p,q= 8  was  determined  by  trial  and  error.  A  simi¬ 
lar  procedure  for  determining  the  stage  one  model  order  was  used  in  Hannan  et  al.  (1986).  The 
SSLS  4  model  parameters  were:  a  =  (l. 58607,— 0.67426,— 0.17672,0.16647), 

6  =  (0. 19765, -0.44775, -0.24888, 0.17300).  The  appearance  of  the  superimposed  SP 
A/=4,  It,  =  1, it2=l  and  3SLS  modeled  impulse  response,  transfer  function  and  phase  response  were 
visually  indistinguishable.  On  the  basis  of  this  evidence,  it  was  thought  reasonable  to  examine  the 
performance  of  the  3SLS  transfer  function  model  with  AR  observation  noise  that  was  similar  to 
the  observation  noise  in  the  BJ  series  J  data. 

The  model  that  we  used  to  synthesize  data  for  the  Monte  Carlo  simulations  is  a  variation  of 
the  model  of  the  BJ  series  J  data.  The  input  data  for  the  simulation  was  the  Box-Jenkins  Series  J 
input  data.  For  the  first  set  of  trials,  the  added  noise  was  a  stochastic  version  of  an  AR4  model  of 
the  residual  noise  from  the  SP  fit  to  the  Box-Jenkins  series  J  data.  The  a,b  coefficients  of  the 
noiseless  simulation  model  were  a  =  (l. 66288,— 0.64256,— 0.30648,0.22377) 

6  =  (-0.83218, -0.47872, -0.24869, 0.12831).  The  AR4  model  coefficients  were 

c  =  (1.69069, -0.69023,-0.28565, 0.022507),  <r2=0.284.  The  (biased)  SP  model  parameters  fitted  to 
a  noiseless  version  of  that  data  were  a  =  (l. 73481, -0.21888, -0.92282, 0.40184), 
fc  =  (-0.03550, -0.05907, 0.04443, 0.05015).  The  vector  of  hyperparameters  was 

r  =  (0.000062,0.000004,0.000027,0.000003).  Such  small  values  should  not  be  surprising  because  to 
within  roundoff  errors,  the  noiseless  data  is  exactly  an  AR4MA4  model. 

Results  of  the  statistical  properties  of  25  replications  from  the  SP  and  SSLS  models  for 
n  =  296,  AR4,  data  points  are  shown  in  Table  2.  The  output  deta  is  regressed  partially  upon  itself 
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so  the  SSLS  procedure  as  well  as  the  SP  procedure,  will  yield  biased  coefficient  estimates.  (The 
magnitude  of  the  biases  is  model  dependent.)  The  bias  of  the  SP  modeled  parameters  is  defined 
as  the  difference  between  the  mean  SP  parameters  and  the  tero  added  noise  SSLS  model  parame¬ 
ters.  The  standard  deviation  and  bias  errors  are  comparable  for  both  the  the  SP  and  SSLS  pro¬ 
cedures.  (The  standard  deviation  of  the  SP  modeled  estimated  b  polynomial  parameters  are  actu¬ 
ally  somewhat  smaller  than  those  for  the  SSLS  procedure.) 


Table  2  AIC’s, 

SP  k  3SLS  M-4  Models 

SP  M  =  4 

SSLS  M  =  4 

param 

mean 

std.  dev. 

bias 

mean 

std.  dev. 

bias 

arl 

-1.6436 

0.0541 

-0.0820 

-1.6641 

0.0554 

-0.0615 

ar2 

-0.7746 

0.109S 

-0.5833 

-0.8S95 

0.1052 

-0.6483 

arS 

-0.1577 

0.1082 

-0.7836 

-0.0911 

0.1054 

0.8502 

ar4 

0.2040 

0.0524 

-0.2029 

0.1775 

0.0538 

-0.2294 

mal 

-0.0404 

0.0217 

-0.0050 

-0.0596 

0.0484 

0.0243 

ma2 

-0.0284 

0.0146 

0.0310 

•0.0181 

0.0893 

0.0414 

maS 

-0.0161 

0.0126 

0.0281 

-0.0078 

0.0861 

0.0519 

ma4 

-0.0192 

0.0106 

-0.0200 

-0.0511 

0.0494 

-0.1019 

Figure  7  is  an  illustration  of  the  mean  and  plus  and  minus  on  sigma  of  results  estimated 
from  the  Monte  Carlo  trials.  The  illustrations  correspond  to  the  simulation  results  reported  in 
Table  2.  From  Figure  7  it  appears  that  the  overall  mean  square  error  in  transfer  function  estimate 
is  slightly  smaller  for  the  SP  than  for  the  SSLS  method.  The  similarity  of  the  SP  and  SSLS  simu¬ 
lation  results  is  compatible  with  the  similarity  of  performance  of  those  models  on  the  BJ  series  J 
data. 


INSERT  FIGURE  7  HERE  SP  AND  SSLS  TRANSFER  FCN  MEAN  AND  STD  DEVS. 

The  order  4a, 6  polynomial  SSLS  model  tends  to  be  overparameterised.  In  order  to  verify 
that  the  SSLS  statistical  results  shown  in  Table  2  were  representative,  an  ARMAX  2,2,2  model 
was  also  simulated  and  modeled  by  the  SSLS  and  SP  M=A  model  procedures.  As  before,  the  input 
was  the  Box-Jenkins  series  J  input  data.  The  additive  stochastically  modeled  noise  was  an  A/?, 
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version  of  the  midul  of  the  SP  modeled  series  J  date.  The  simulation  model  parameters  were 
a*(1.46778, -0.55 192), 8*(-0.03586, -0.06788),*  *(1.58778, -0.72461),  <7**0.078.  The  statistical 
results  of  25  replications  of  fitting  the  SP  and  SSLS  models  to  the  trial  data  are  shown  in  Table  S. 
Here,  the  bias  of  the  SP  modeled  parameters  is  defined  as  the  difference  between  the  mean  and  sero 
added  noise  SP  model  parameters.  The  standard  deviation  and  bias  of  the  SP  A#»4  model  param¬ 
eters  are  comparable  to  those  in  Table  2.  The  standard  deviation  and  bias  of  the  SSLS 
model  parameters  are  similar  to  those  observed  in  Table  2  for  the  Af*4  model. 


Table  3.  AIC  SP  M«4,  k  3! 

SP  M  =  4 

5LS  M=2  Models 

SSLS  M  =  2 

par  am 

mean 

std.  dev. 

bias 

mean 

std.  dev. 

bias 

arl 

1.5846 

0.0566 

-0.1502 

-1.1562 

0.0386 

-0.1633 

ar2 

-0.7048 

0.1072 

-0.4910 

-0.0095 

0.0356 

-0.1171 

arS 

-0.0323 

0.0984 

0.8905 

- 

- 

* 

art 

0.0486 

0.0415 

-0.3532 

- 

* 

* 

mal 

-0.0443 

0.0222 

-0.0088 

-0.0095 

0.0406 

0.0264 

ma2 

•0.0362 

0.03 83 

0.0229 

-0.0918 

0.0436 

-0.0231 

maS 

-0.0001 

0.0243 

•0.0451 

- 

* 

• 

mat 

-0.0299 

0.0236 

0.0600 

- 

- 

- 

We  recall  that  the  added  noise  for  the  Table  2  data  was  ARt  and  that  for  Table  S  was  ARt- 
The  consistency  of  the  tabulated  results  for  the  SP  model  in  Tables  2  and  3  suggest  that  the  SP 
model  is  reasonably  robust  with  respect  to  noise  color  and  noise  model  order. 

Finally,  we  show  results  of  simulation  studies  with  added  MA  observation  noise.  The  pri¬ 
mary  objects  of  interests  in  these  computational  experiments  were  a  comparison  of  the  SP  and 
SSLS  modeling  performance  and  the  sensitivity  of  the  SP  transfer  function  modeling  method  to 
choices  of  model  order  M.  the  differential  orders  kjk2,  the  observation  noise  level  and  the  sensi¬ 
tivity  to  data  length. 

The  model  for  these  simulations  was  a  slight  variation  of  the  model  for  the  AR  observation 
noise  simulations.  The  superimposed  impulse  response,  transfer  function  and  phase  function  of  the 
SSLS  and  SP  M=4,M=5  and  M*  10  on  the  noise  free  data  were  visually  indistinguishable.  In  the 
first  stage  of  the  SSLS  method  a  model  order  of  20  was  used.  The  stochastic  input  signal  was  an 
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AR4  version  of  tke  iapat  of  the  series  J  data.  The  SSLS  W»4  noise  free  data  a,  6  polynomials 
wnpe:  •  -  {1.72566,-0.19127,-0.94126,0.40697},  b  -  {-0  03550,-0.05940,0.04412,0.05070}. 

The  observation  noise  was  an  MAt  model,  (inadvertently  MAS,  instead  of  MAt).  The  noise 
■oM  parameters  were  e  -  {1.0, 0.6, 0.3, 0.1}.  This  model  has  a  fairly  flat  psd  spectrum.  A  sample 
run  of  the  MAl  noise  yielded  an  ARX 0  AIC-AR  model  of  that  noise. 

Monte  Carlo  trials  with  n “296  and  »»T8  were  done.  Superimposed  output  and  output  plus 
noise  and  a  displaced  version  of  the  noise  are  shown  in  Figure  8  for  typical  sample  trials  of  data 
lengths  » -296  and  it -78.  Figure  8  illustrates  the  noise  level  of  the  experiments  and  how  dramati¬ 
cally  short  the  n-78  data  span  is. 

Differentia)  orders  kx  and  best  SP  A/=5,Af  =  10  and  M=  20  transfer  function  models  were 
ftted  to  a  sample  trial  of  the  n-296  data  with  the  following  results: 
M-5,  *,-1,  *,-l,  A/C=- 1567.545,  A/-10,  *,=4,  *s=l,  A/C—1609.348,  Af-20,  *,=4,  *,-3, 
A/C- -1557.410.  The  similarity  of  the  AIC  values  suggest  that  these  modeb  will  have  similar 
properties. 

The  computed  results  of  the  MAj  noise  simulation  runs  for  the  SSLS  M-4, 
SP  M-10,  fc,-l,*J=4  and  SP  Af=  20,  *,=2,*j=4,  SP  M= 20,  *,  =  1,^=3,  SP  Af=20,  *,-S,*j=5  and 
an  SP  A/=20,  i,  =  2,ij=4  model  of  data  with  observation  noise  variance  9  times  that  of  the  previ¬ 
ous  three  simulation  trials  are  shown  in  Figure  9A,B,C,D,E.F.  Tabulated  values  of  the  mean  and 
standard  deviation  of  the  parameter  estimates  for  the  SSLS  and  SP  Af=10,  Af=20  1,  =2,ij=4 
simulation  trials  are  shown  in  Table  4.  The  SSLS  and  SP  Af=10  results  transfer  function  and 
phase  function  results  are  reasonable  similar.  The  relative  jaggedness  of  the  impulse  response  for 
the  SSLS  model  indicates  that  a  lower  order  SSLS  mode)  might  be  more  suitable  for  this  data. 
The  SP  A/=20,  kl,kJ=4  model  results  appear  very  similar  to  the  results  for  the  M=  10  modeb.  The 
results  of  the  Monte  Carlo  runs,  with  SP  modeb  Af-20,l:,-l,fcj=S  and  M=20,t,=S,fcj=5,  (Figures 
9D,E),  are  very  similar  to  those  shown  for  the  A/=20,i,=2,ij=4,  (Figure  9C),  results.  Those  illus¬ 
trations  show  the  relative  insenstivity  of  SP  modeling  to  variations  in  the  differential  orders  At ,lr2. 
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As  expected,  the  results  in  Figure  9F  for  the  SP  M=20  mode)  the  9  times  variance  trials  show 
more  dispersion  than  the  results  for  the  Figure  9C  trials.  The  graphical  results  in  Figure9  are  com¬ 
patible  with  the  performance  of  the  SSLS  and  SP  models  of  the  Box-Jenkins  series  J  data. 


Table  4. 

AIC,  3SLS  M=4  &  SP  M=  10,20,20  Models 

SSLS  M=4 

SP  M 

=  10 

SP  M 

=20 

SP  M=20,  VAR  x  9 

■a 

s.d. 

mean 

s.d. 

mean 

s.d. 

mean 

s.d. 

arl 

0.7217 

0.0469 

0.6420 

0.0356 

0.6100 

H 

ar2 

0.0001 

0.0588 

-0.0109 

■  jXlVV£jH| 

-0.0154 

1 

WITWTWWSt 

arS 

0.0630 

-0.0013 

y  1 

-0.0000 

1 

0.0068 

ar4 

-0.0041 

0.0397 

0.0031 

1 

0.0002 

0.0001 

-0.0001 

0.0005 

ar5 

- 

- 

0.0016 

0.0001 

0.0001 

0.0001 

0.0002 

mat 

0.0235 

-0.0363 

-0.0437 

0.0248 

-0.0547 

0.0463 

ma2 

0.0498 

-0.0958 

0.0366 

-0.0906 

0.0366 

-0.0907 

0.0621 

maS 

0.0060 

0.0557 

-0.0798 

0.0318 

-0.0884 

0.0308 

-0.0675 

0.1091 

ma4 

-0.1526 

0.0333 

-0.0719 

0.0277 

-0.0674 

0.0252 

-0.0879 

0.0608 

ma5 

- 

- 

-0.0467 

0.0243 

-0.0628 

0.0253 

-0.0704 

0.0539 

In  Table  4  only  the  results  for  the  first  5  a,b  polynomial  coefficients  are  shown.  That  is 
sufficient  to  understand  the  results  of  the  Monte  Carlo  trials.  The  standard  deviations  for  the 
3SLS  and  M=  10  SP  trials  are  comparable.  Those  results  are  also  comparable  to  those  shown  in 
Table  2  for  the  AR  noise  trials.  The  a  polynomial  coefficients  for  the  SP  trials  tend  to  taper 
quickly  and  correspondingly,  the  standard  deviations  for  those  coefficients  tend  to  be  smaller  than 
the  standard  deviations  for  the  larger  valued  coefficients.  The  values  of  the  b  polynomial 
coefficients  do  not  taper  and  correspondingly  neither  do  the  standard  deviations  of  those  coefficient 
estimates. 


The  final  Monte  Carlo  trials  arc  for  the  shorter  length  n=78  data.  An  exploratory  computa¬ 
tion  on  a  single  trial  of  an  iV=78  data  SP  A#*  10  model,  (that  yields  an  N= 68  data  points  for 
actual  modeling),  indicated  that  the  A,  =  2,4c2=  1  model  was  optimum  for  that  data.  No  experiment 
was  performed  for  the  SP  M=5  model.  Graphical  results  for  25  SP  M=5,k1  =  l,kJ=l  and  SP 
Af  =—  10,Aj  “  1, A2=- 1  models  MA  noise  Monte  Carlo  trials  are  shown  in  Figure  10.  Numerical  results 
corresponding  to  those  shown  in  Figure  10  are  in  Table  5.  (Only  results  for  the  first  6  polynomial 
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coefficients  are  show*.)  The  standard  deviations  for  the  n-78  data  coefficients  are  in  general 
larger  than  those  for  the  *=296  data  coefficients.  The  SP  M- 10  *-78  data  parameter  mean 
valves  and  standard  deviations  indicate  the  same  kind  of  tapering  effects  as  was  observed  for  the 
n-296  data.  The  graphical  resalts  for  the  SP  transfer  function  modeled  n=296  and  n«78  data 
look  very  similar.  This  is  additional  evidence  to  support  the  conjecture  that  Bayesian  smoothness 
priors  transfer  function  modeling  can  yield  reliable  transfer  function  estimation  with  relatively 
short  duration  data. 


Table  5.  AIC,  3SLS  M~4  k  SP  M- 10,20,20  Models 

SP  M  =  4 

SP  M  =  10 

par  am 

mean 

std-  dev. 

mean 

std.  dev. 

arl 

0.7217 

0.0848 

0.6623 

0.0816 

ar2 

0.0196 
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0.0537 

0.0671 
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0.0382 

0.0678 

0.0186 

0.0605 
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0.0164 

0.0318 
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-0.0151 

0.0297 

0.0153 
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- 

- 

0.0062 
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-0.0416 

0.0551 
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-0.0772 

0.0344 

-0.0971 

0.0753 

maS 

-0.0735 

0.0335 

-0.0663 

0.0426 
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-0.0521 

0.0322 

-0.0616 

0.0352 

maS 

-0.0512 

0.0342 

-0.0330 

0.0382 

ma6 

- 

- 

-0.0125 

0.0162 
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4.  INTERPRETATION  OP  RESULTS,  SUMMARY  AND  DISCUSSION 

A  smoothness  priors  approach  to  the  problem  of  transfer  function  estimation  was  demon¬ 
strated.  The  smoothness  priors  are  stochastic  constraints  on  the  parameters  of  the  linear  model. 
The  Bayesian  smoothness  priors  procedure  is  one  possible  way  of  utilising  the  information  in  the 
likelihood  function.  The  critical  computation  in  the  Bayesian  smoothness  priors  approach  is  the 
likelihood  of  the  Bayesian  model.  The  likelihood  is  used  as  an  objective  measure  of  the  goodness  of 
a  model  and  the  hyperparameters  which  maximise  the  likelihood  are  determined  by  a  gradient 
search  method. 

Some  of  the  consequences  of  the  smoothness  priors  approach  to  transfer  function  estimation 

are: 

1)  The  complex  multiple  model  orders  selection  procedures  of  other  transfer  function  estimation 
methods  are  obviated  by  the  smoothness  priors  method.  Instead,  the  specification  of  the  values  of 
model  order,  M,  and  the  differential  orders,  kltkt,  appear  not  to  be  very  critical  in  determining  the 
impulse  response,  transfer  function  and  phase  properties  of  the  model. 

2)  Least  squares  and  conventional  maximum  likelihood  are  abandoned  as  a  criterion  for  modeling 
data.  Instead,  constrained  least  squares  criterion  or  equivalently  penalised  likelihood  methods  are 
used  in  the  smoothness  prior  method. 

3)  The  smoothness  priors  model  will  in  general  not  be  as  parsimonious,  in  the  number  of  model 
parameters,  as  the  conventional  least  squares  or  maximum  likelihood  transfer  function  estimation 
methods. 

4)  Smoothness  priors  modeling  will  tend  to  be  more  robust  with  respect  to  the  specification  of  the 
observation  noise  color  than  the  conventional  methods. 

5)  Smoothness  priors  modeling  will  tend  to  yield  reliable  parameter  estimates  with  shorter  length 
data  than  conventional  methods. 

The  relatively  large  model  order  non-parsimoniousness  property  of  the  Bayesian  smoothness 
priors  model  can  be  justified.  In  any  alternative  Bayesian  transfer  function  modeling,  the  Bayesian 


model  would  be  an  average  of  the  transfer  function  computed  by  models  of  the  different  a,b,c  poly* 
nomial  model  orders  weighted  by  the  likelihood  and  the  priors  on  model  orders.  The  Bayesian 
transfer  function  model  would  have  contributions  from  the  largest  a,b,c  polynomial  model  orders 
considered.  As  a  result,  the  overall  Bayesian  model  orders  would  be  the  identical  to  the  orders  of 
that  largest  orders  model.  Thus,  the  Bayesian  procedure  will  have  a,b,c  polynomial  model  orders 
larger  than  the  model  selected  say  by  the  minimum  AIC  procedure.  By  the  specification  of  a  large 
model  order  and  smoothness  priors  constraints  on  the  model  parameters,  the  smoothness  priors 
method  achieves  the  effect  of  the  alternative  Bayesian  methods  of  averaging  the  computational 
results  of  different  order  models. 

The  use  of  frequency  domain  priors  is  a  relatively  new  idea.  The  class  of  frequency  domain 
priors  that  we  used  lack  a  definite  physical  interpretation.  They  seem  reasonable  because  they 
penalise  the  higher  order  polynomial  coefficients  with  increasing  weights.  As  expected,  the  SP 
model  properties  are  increasingly  smooth  with  increasing  kth  derivative  constraints.  Also,  the 
overall  optimum  solution  is  not  necessarily  the  smoothest  solution.  An  important  property  of  our 
frequency  domain  priors  is  that  they  permit  the  problem  of  transfer  function  estimation  to  be  cast 
within  the  framework  of  the  linear  model. 

The  Monte  Carlo  results  suggest  that  the  smoothness  priors  method  of  transfer  function  esti¬ 
mation  achieves  comparable  statistical  performance  with  the  asymptotically  efficient  procedures. 

In  summary,  the  smoothness  priors  method  of  transfer  function  estimation  appears  to  yield 
statistical  performance  results  that  are  comparable  and  perhaps  superior  to  other  well  known 
methods  while  enjoying  a  more  forgiving  parameter  computational  search  procedure  than  the 
search  procedures  of  other  methods.  In  general,  the  results  shown  attest  to  the  flexibility  of  the 
Bayesian  modeling  approach. 
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LEGENDS 


FIGURE  1.  Box-Jenkins  Series  J  Gas  Furnace  Data 

Figure  2.  Original  and  Modeled  Data  Superimposed  and  Difference  Data  Smoothness  Priors 
Model,  n-296 

FIGURE  3.  Impulse  Response,  Amplitude  and  Phase  Response  and  Residual  PSD  versus  Fre¬ 
quency,  Superimposed  Box-Jenkins  and  Smoothness  Priors  Model  Results,  (Dotted  Lines). 

FIGURE  4.  Superimposed  Smoothness  Priors  Impulse  Response,  Amplitude  Response,  Phase 
Response  and  Noise  Spectrum  for 

A:  n-296  and  n=200  data  points,  (dotted  lines).  B:  n=296  and  n=75  data  points,  (dotted  lines). 

FIGURE  5.  Superimposed  Impulse  Response,  Amplitude  Response,  Phase  Response  and  Noise 
Spectrum: 

A:  SP  Af=>4  and  M=10  models.  B:  SP  M=4  and  M-20  models. 

C:  SP  M= 4  and  A/=30  models.  D:  SP  Af=4  and  LS  M=  5  models. 

E:  SP  Af=4  and  LS  A/-IO  models.  F:  SP  4  and  LS  Af=20  models. 

FIGURE  6.  Superimposed  Smoothness  Priors  Impulse  Response,  Amplitude  Response,  Phase 
Response  and  Noise  Spectrum  for  SP  Af=4  kt  =  l,k ?=I  and  kl=9,kJ=9  (dotted  lines)  models. 

FIGURE  7.  Mean  and  Plus  and  Minus  One  Standard  Deviation  Impulse  Response,  Amplitude 
Response,  Phase  Response.  AR  Observation  Noise: 

A:  Smoothness  Priors  Model.  B:  SSLS  Model. 

FIGURE  8.  Input.  Output  and  Ouput  Plus  Noise: 

A:  N* 296,  Input,  Noise,  Superimposed  output  and  output  plus  noise,  Superimposed  output  and 
output  plus  larger  variance  noise,  (in  descending  order). 

B:  N-78,  Superimposed  output  and  output  plus  noise.  Noise. 


FIGURE  9.  Mean  and  Plus  and  Minus  One  Standard  Deviation  Impulse  Response.  Amplitude 
Response  and  Phase  Response.  MA  Observation  Noise.  A’  =  296: 

A:  SSLS  A/=  4  Model.  B:  SP  M=10,  *,  =  2,*,  =  4  Model. 

C:  SP  M= 20.  L,  =  l,Jb2=S  Model.  D:  SP  M=20.  L,  =  3,fcj=5  Model. 

E:  SP  A/=20,  ki  =  2,kI=4  Model.  F:  SP  M=20,  kt=2,k2=4  Model,  Variance  x  9. 

FIGURE  10.  Mean  and  Pius  and  Minus  One  Standard  Deviation  Impulse  Response,  Amplitude 
Response  and  Phase  Response,  MA  Observation  Noise,  N  =  78: 

A:  SP  M=  5  *,=2, Jfc2=l  Model.  B:  SP  M=10 ,  i,  =  2,*2=l  Model. 
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